# File Name:  RegsMonteCarloStudy_jpp.R
# Data:       MCRegsData.RData
# Purpose:    Create a simulated data set based on the Shor-McCarty Data to estimate regressions when the lawmaking model is known. 
#             Simulation fills in missing data with the data from the previous state-period.
# Output:     MCreg.RData; MCregs_three_fifths.RData; MCregs_unanimity.RData; MCsummary_fd; MCsummary_level;


setwd("C:/Users/MDR/Dropbox/Vanderbilt/Fall 2013/Project/Charter/ClintonRichardson Draft/MonteCarlo")

load("MCRegsData.RData")

library("doParallel")
library("foreach")

state<-unique(MCData$state) #get vector of state names for the sim function

#get a list of missing state-years
for(i in 1:nrow(MCData)){
  if(is.na(MCData$hou_dem.avg[i])==T){
    print(MCData$state.year[i])
  }
}

#write a function to replace the missing data
replace<- function(missing,replace){
  
  MCData$hou_dem.avg[MCData$state.year==missing]<-MCData$hou_dem.avg[MCData$state.year==replace]
  MCData$hou_dem.sd[MCData$state.year==missing]<-MCData$hou_dem.sd[MCData$state.year==replace]
  
  MCData$hou_rep.avg[MCData$state.year==missing]<-MCData$hou_rep.avg[MCData$state.year==replace]
  MCData$hou_rep.sd[MCData$state.year==missing]<-MCData$hou_rep.sd[MCData$state.year==replace]
  
  MCData$rep.size[MCData$state.year==missing]<-MCData$rep.size[MCData$state.year==replace]
  return(MCData)
}
###########
#recode missing state-years
###########
MCData<-replace(missing="AK2010", replace="AK2008")
MCData<-replace(missing="AR2010", replace="AR2008")
MCData<-replace(missing="CT2010", replace="CT2008")
MCData<-replace(missing="DE2010", replace="DE2008")
MCData<-replace(missing="FL2010", replace="FL2008")
MCData<-replace(missing="FL2013", replace="FL2012")
MCData<-replace(missing="GA2010", replace="GA2008")
MCData<-replace(missing="GA2013", replace="GA2012")
MCData<-replace(missing="HI2010", replace="HI2008")
MCData<-replace(missing="IA2013", replace="IA2012")
MCData<-replace(missing="IN2010", replace="IN2008")
MCData<-replace(missing="KY1999", replace="KY1998")
MCData<-replace(missing="KY2010", replace="KY2008")
MCData<-replace(missing="KY2012", replace="KY2008")
MCData<-replace(missing="KY2013", replace="KY2008")
MCData<-replace(missing="MA2010", replace="MA2008")
MCData<-replace(missing="MA2012", replace="MA2008")
MCData<-replace(missing="ME2010", replace="ME2008")
MCData<-replace(missing="MI2010", replace="MI2008")
MCData<-replace(missing="MN2013", replace="MN2012")
MCData<-replace(missing="MO2010", replace="MO2008")
MCData<-replace(missing="MO2012", replace="MO2008")
MCData<-replace(missing="MO2013", replace="MO2008")
MCData<-replace(missing="MT2010", replace="MT2008")
MCData<-replace(missing="ND2010", replace="ND2008")
MCData<-replace(missing="ND2013", replace="ND2012")
MCData<-replace(missing="NM2010", replace="NM2008")
MCData<-replace(missing="OH2010", replace="OH2008")
MCData<-replace(missing="OK2010", replace="OK2008")
MCData<-replace(missing="TN2010", replace="TN2008")
MCData<-replace(missing="TX2013", replace="TX2012")
MCData<-replace(missing="UT2010", replace="UT2008")
MCData<-replace(missing="WV2010", replace="WV2008")

#confirm no missing remain
for(i in 1:nrow(MCData)){
  if(is.na(MCData$hou_dem.avg[i])==T){
    print(MCData$state.year[i])
  }
}

#recode missing MN years for Jessie Ventura - make him a Republican
table(MCData$state,is.na(MCData$gov.rep))

MCData$gov.rep[MCData$state=="MN"&is.na(MCData$gov.rep)==T]<-1

############################
#create data sets for simulating two additional filibuster assumptions
############################

#For assuming each state requires a 3/5 supermajority to invoke cloture


############################
#check for filibuster pivots more extreme than the veto pivot in MCData
############################

for(i in 1:49){
  if(unique(MCData$buster_norm[MCData$state==state[i]]==1)){
    print(state[i])
    print(c("vp:",unique(MCData$vp[MCData$state==state[i]])))
    print(c("fp:",unique(MCData$fp_norm[MCData$state==state[i]])))
  }
}

#Vermont is the only case where fp>vp
table(MCData$vp[MCData$state=="VT"])
MCData$vp[MCData$state=="VT"]<-0.75
table(MCData$vp[MCData$state=="VT"])

############################
#Create data set for simulation assuming all states have 3/5 supermajority requirement to invoke cloture
############################
MC35<-MCData 

#check for filibuster pivots more extreme than the veto pivot
table(MC35$state[MC35$vp<3/5]) #affects 6 states
MC35$vp[MC35$vp<3/5]<-3/5

############################
#Create data set for simulation assuming certain states require unanimity
############################

MCunam<-MCData

#code veto pivot percentile as 1 in states with no way to shut off debate so that unanimity is required to change policy
MCunam$vp[MCunam$fp==1]<-1 

#Note: Vermont upper chamber is coded as having no way to shut off debate
#and Vermont lower chamber requires 3/4 to move the previous question.
#Therefore, it will have 3/4 fp in model in main text and requires unanimity if no norm of limited debate is assumed.

#check for filibuster pivots more extreme than the veto pivot
for(i in 1:49){
  if(unique(MCunam$buster[MCunam$state==state[i]]==1)){
    print(state[i])
    print(c("vp:",unique(MCunam$vp[MCunam$state==state[i]])))
    print(c("fp:",unique(MCunam$fp[MCunam$state==state[i]])))
  }
} #No recoding needed

###############
#write a function to model each session for each theory
#function takes a vector ideal point averages and a vector standard deviations of ideal points for each party,
# a vector of party indicators for the governor (republican =1, democrat=0), a vector of the size of the republican caucus, 
# a starting status quos, and a number of sessions to model.
#takes a veto pivot (vp) and a filibuster pivot (fp)
#Other vectors should be of length=sessions

lawmaking<-function(dem.avg,dem.sd,rep.avg,rep.sd,gov.rep,length.rep,vp,fp,sq,sessions){

  #create vectors of preferences
  xmaj<-vector(length=sessions)
  xmedian<-vector(length=sessions)
  maj.name<-vector(length=sessions)

  #session counter
  session<-vector(length=sessions)
  
  #create vectors of policy outcomes
  p.pivot<-vector(length=sessions)
  p.pivot.small<-vector(length=sessions)
  p.median<-vector(length=sessions)
  p.gatekeep<-vector(length=sessions)
  p.positive<-vector(length=sessions)
  p.exec<-vector(length=sessions)

  #set the status quos
  sq.pivot<- sq
  sq.pivot.small<-sq
  sq.median<-sq
  sq.gatekeep<-sq
  sq.positive<-sq
  sq.exec<-sq

  for(i in 1:sessions) {
    n.rep<-length.rep[i] #get the number of republicans for that state-year
    n.dem<-100-n.rep

    dems<-rnorm(n.dem,mean=dem.avg[i],sd=dem.sd[i]) #draw the dems
    reps<-rnorm(n.rep,mean=rep.avg[i],sd=rep.sd[i]) #draw the reps
    
    #set the majority party median
    if(n.dem > n.rep) maj<-quantile(dems,.5) else maj<-quantile(reps,.5)
    if(n.dem > n.rep) maj.name[i]<-"dems" else maj.name[i]<-"reps"
    xmaj[i]<-maj
    
    #set the chamber median
    all<-c(dems,reps)
    med<-quantile(all,.5)
    xmedian[i]<-quantile(all,.5)
    
    #set the pivots - based on party of the governor and state's veto provision;
    if(gov.rep[i]==1){
      fr<-quantile(all,vp)
      fl<-quantile(all,(1-fp))
    }
    
    if(gov.rep[i]==0){
      fr<-quantile(all,fp)
      fl<-quantile(all,(1-vp))
    }
    
    #pivot model
    if(sq.pivot <= 2*fl - med)                    p.pivot[i]<- med
    if(sq.pivot > (2*fl - med) & sq.pivot < fl)   p.pivot[i]<- (2*fl-sq.pivot)
    if(sq.pivot >= fl & sq.pivot <= fr)          p.pivot[i]<- sq.pivot
    if(sq.pivot <= (2*fr - med) & sq.pivot > fr)  p.pivot[i]<- (2*fr-sq.pivot)
    if(sq.pivot > (2*fr - med))                   p.pivot[i]<- med
    
    #pivot model - no filibuster
    if(gov.rep[i]==1){
      if(sq.pivot.small <= med)                                 p.pivot.small[i]<- med
      if(sq.pivot.small > med & sq.pivot.small <= fr)           p.pivot.small[i]<- sq.pivot.small
      if(sq.pivot.small <= (2*fr - med) & sq.pivot.small > fr)  p.pivot.small[i]<- (2*fr-sq.pivot.small)
      if(sq.pivot.small > (2*fr - med))                         p.pivot.small[i]<- med
    }
  
    if(gov.rep[i]==0){
      if(sq.pivot.small >= med)                                 p.pivot.small[i]<- med
      if(sq.pivot.small < med & sq.pivot.small >= fl)           p.pivot.small[i]<- sq.pivot.small
      if(sq.pivot.small >= (2*fl - med) & sq.pivot.small < fl)  p.pivot.small[i]<- (2*fl-sq.pivot.small)
      if(sq.pivot.small < (2*fl - med))                         p.pivot.small[i]<- med
    }
  
    #median voter
    p.median[i] <- med
    
    #negative agenda power
    if(maj > med){  
      if(sq.gatekeep < med)                                 p.gatekeep[i]<- med
      if(sq.gatekeep >= med & sq.gatekeep <= (2*maj - med)) p.gatekeep[i]<- sq.gatekeep
      if(sq.gatekeep > (2*maj - med))                       p.gatekeep[i]<- med
      }
    
    if(maj < med){  
      if(sq.gatekeep < (2*maj - med))                       p.gatekeep[i]<- med
      if(sq.gatekeep >= (2*maj - med) & sq.gatekeep <= med) p.gatekeep[i]<- sq.gatekeep
      if(sq.gatekeep > med)                                 p.gatekeep[i]<- med
      }
    
    #positive agenda power
    if(maj > med){ 
      if(sq.positive <= (2*med-maj))                    p.positive[i]<- maj
      if(sq.positive > (2*med-maj) & sq.positive < med) p.positive[i]<- (2*med-sq.positive)
      if(sq.positive >= med & sq.positive <= maj)       p.positive[i]<- sq.positive
      if(sq.positive > maj)                             p.positive[i]<- maj
    }
    
    if(maj < med){ 
      if(sq.positive < maj)                               p.positive[i]<- maj
      if(sq.positive >= maj & sq.positive <= med )        p.positive[i]<- sq.positive
      if(sq.positive > med & sq.positive <= (2*med-maj))  p.positive[i]<- (2*med-sq.positive)
      if(sq.positive > (2*med-maj))                       p.positive[i]<- maj
      }
    
    #executive agenda setter; assume governor's ideal point is the party median;
    
    #set the minority party median for executive agenda setting model
    if(n.dem > n.rep) min<-quantile(reps,.5) else min<-quantile(dems,.5)
    
    #set executive ideal point to majority party median in united government
    if( (gov.rep[i]==1 & maj.name[i]=="reps")|(gov.rep[i]==0 & maj.name[i]=="dems") ) exec<- maj
    
    #set executive ideal point to minority party median in divided government
    if( (gov.rep[i]==0 & maj.name[i]=="reps")|(gov.rep[i]==1 & maj.name[i]=="dems") ) exec<- min
      
    if(exec > med){ 
      if(sq.exec <= (2*med-exec))                  p.exec[i]<- exec
      if(sq.exec > (2*med-exec) & sq.exec < med)   p.exec[i]<- (2*med-sq.exec)
      if(sq.exec >= med & sq.exec <= exec)         p.exec[i]<- sq.exec
      if(sq.exec > exec)                           p.exec[i]<- exec
    }
      
    if(exec < med){ 
      if(sq.exec < exec)                           p.exec[i]<- exec
      if(sq.exec >= exec & sq.exec <= med )        p.exec[i]<- sq.exec
      if(sq.exec > med & sq.exec <= (2*med-exec))  p.exec[i]<- (2*med-sq.exec)
      if(sq.exec > (2*med-exec))                   p.exec[i]<- exec
    }

    #set status quo to policy outcome from current session
    sq.pivot<- p.pivot[i]
    sq.pivot.small<- p.pivot.small[i]
    sq.median<- p.median[i]
    sq.gatekeep<- p.gatekeep[i]
    sq.positive<- p.positive[i]
    sq.exec<- p.exec[i]
  }
  
  #assemble the data
  data<-cbind(xmaj,xmedian,p.pivot,p.pivot.small,p.median,p.gatekeep,p.positive,p.exec,gov.rep)
  
  return(data) 
}


########################################################
#create a function to assemble a simulated data set using for a vector of state abbreviations and a starting status quo;

sim<-function(state){
  
  #create a matrix for the simulated data
  data<-matrix(NA,nrow=1,ncol=11)
  colnames(data)<-c("xmaj","xmedian","p.pivot","p.pivot.small","p.median","p.gatekeep","p.positive","p.exec","gov.rep","session","st")
  
  #Loop to model each state session and create a simulated data set
  for(j in 1:length(state)){
    
    sq<-runif(1,min=-3,max=3)
    
    #Get data for the state being modeled
    dem.avg<-MCData$hou_dem.avg[MCData$state==state[j]]
    dem.sd<-MCData$hou_dem.sd[MCData$state==state[j]]
    rep.avg<-MCData$hou_rep.avg[MCData$state==state[j]]
    rep.sd<-MCData$hou_rep.sd[MCData$state==state[j]]
    length.rep<-MCData$rep.size[MCData$state==state[j]]
    gov.rep<-MCData$gov.rep[MCData$state==state[j]]
    vp<-unique(MCData$vp[MCData$state==state[j]])
    fp<-unique(MCData$fp_norm[MCData$state==state[j]])
    sessions<-length(dem.avg)
    
    #get the simulated sessions for the state
    x<-lawmaking(dem.avg=dem.avg,dem.sd=dem.sd,rep.avg=rep.avg,rep.sd=rep.sd,gov.rep=gov.rep,length.rep=length.rep,vp=vp,fp=fp,sq=sq,sessions=sessions)
   
    st<-rep(state[j],length=sessions)#create a state label
    
    session<-MCData$session[MCData$state==state[j]] #get session number for regressions
    
    x<-cbind(x,session,st) #combine the state data, state label, and session number 
    
    data<-rbind(data,x) #add the state to the full data set
  }
  
  data<-as.data.frame(data,stringsAsFactors=F) #convert to a data frame
  
  #convert variables to numeric
  for(t in 1:(ncol(data)-1)){
    data[,t]<-as.numeric(data[,t])
  }
  data<-data[-1,] #remove the first row of NA's used to start the data 
  return(data)
}

##########################################
#Write a function to run regressions and store coefficients and p-values

level.reg<-function(data){

names<-c("pivot.med.coef","pivot.med.p","pivot.maj.coef","pivot.maj.p","pivot.gov.coef","pivot.gov.p",
         "pivot.sm.med.coef","pivot.sm.med.p","pivot.sm.maj.coef","pivot.sm.maj.p","pivot.sm.gov.coef","pivot.sm.gov.p",
         "median.med.coef","median.med.p","median.maj.coef","median.maj.p","median.gov.coef","median.gov.p",
         "gatekeep.med.coef","gatekeep.med.p","gatekeep.maj.coef","gatekeep.maj.p","gatekeep.gov.coef","gatekeep.gov.p",
         "positive.med.coef","positive.med.p","positive.maj.coef","positive.maj.p","positive.gov.coef","positive.gov.p",
         "exec.med.coef","exec.med.p","exec.maj.coef","exec.maj.p","exec.gov.coef","exec.gov.p")
  
#create a dataframe to store regression coefficients
mc<-matrix(NA,nrow=1,ncol=length(names))
colnames(mc)<-names
  
#estimate models
pivot<-lm(p.pivot~xmedian+xmaj+gov.rep+session+as.factor(st)-1,data=data)
pivot.small<-lm(p.pivot.small~xmedian+xmaj+gov.rep+session+as.factor(st)-1,data=data)
median<-lm(p.median~xmedian+xmaj+gov.rep+session+as.factor(st)-1,data=data)
gatekeep<-lm(p.gatekeep~xmedian+xmaj+gov.rep+session+as.factor(st)-1,data=data)
positive<-lm(p.positive~xmedian+xmaj+gov.rep+session+as.factor(st)-1,data=data)
exec<-lm(p.exec~xmedian+xmaj+gov.rep+session+as.factor(st)-1,data=data)

#store coefficients and p-values
mc[1,"pivot.med.coef"]<-summary(pivot)$coefficients["xmedian",1]
mc[1,"pivot.med.p"]<-summary(pivot)$coefficients["xmedian",4]
mc[1,"pivot.maj.coef"]<-summary(pivot)$coefficients["xmaj",1]
mc[1,"pivot.maj.p"]<-summary(pivot)$coefficients["xmaj",4]
mc[1,"pivot.gov.coef"]<-summary(pivot)$coefficients["gov.rep",1]
mc[1,"pivot.gov.p"]<-summary(pivot)$coefficients["gov.rep",4]


mc[1,"pivot.sm.med.coef"]<-summary(pivot.small)$coefficients["xmedian",1]
mc[1,"pivot.sm.med.p"]<-summary(pivot.small)$coefficients["xmedian",4]
mc[1,"pivot.sm.maj.coef"]<-summary(pivot.small)$coefficients["xmaj",1]
mc[1,"pivot.sm.maj.p"]<-summary(pivot.small)$coefficients["xmaj",4]
mc[1,"pivot.sm.gov.coef"]<-summary(pivot.small)$coefficients["gov.rep",1]
mc[1,"pivot.sm.gov.p"]<-summary(pivot.small)$coefficients["gov.rep",4]

mc[1,"median.med.coef"]<-summary(median)$coefficients["xmedian",1]
mc[1,"median.med.p"]<-summary(median)$coefficients["xmedian",4]
mc[1,"median.maj.coef"]<-summary(median)$coefficients["xmaj",1]
mc[1,"median.maj.p"]<-summary(median)$coefficients["xmaj",4]
mc[1,"median.gov.coef"]<-summary(median)$coefficients["gov.rep",1]
mc[1,"median.gov.p"]<-summary(median)$coefficients["gov.rep",4]

mc[1,"gatekeep.med.coef"]<-summary(gatekeep)$coefficients["xmedian",1]
mc[1,"gatekeep.med.p"]<-summary(gatekeep)$coefficients["xmedian",4]
mc[1,"gatekeep.maj.coef"]<-summary(gatekeep)$coefficients["xmaj",1]
mc[1,"gatekeep.maj.p"]<-summary(gatekeep)$coefficients["xmaj",4]
mc[1,"gatekeep.gov.coef"]<-summary(gatekeep)$coefficients["gov.rep",1]
mc[1,"gatekeep.gov.p"]<-summary(gatekeep)$coefficients["gov.rep",4]

mc[1,"positive.med.coef"]<-summary(positive)$coefficients["xmedian",1]
mc[1,"positive.med.p"]<-summary(positive)$coefficients["xmedian",4]
mc[1,"positive.maj.coef"]<-summary(positive)$coefficients["xmaj",1]
mc[1,"positive.maj.p"]<-summary(positive)$coefficients["xmaj",4]
mc[1,"positive.gov.coef"]<-summary(positive)$coefficients["gov.rep",1]
mc[1,"positive.gov.p"]<-summary(positive)$coefficients["gov.rep",4]

mc[1,"exec.med.coef"]<-summary(exec)$coefficients["xmedian",1]
mc[1,"exec.med.p"]<-summary(exec)$coefficients["xmedian",4]
mc[1,"exec.maj.coef"]<-summary(exec)$coefficients["xmaj",1]
mc[1,"exec.maj.p"]<-summary(exec)$coefficients["xmaj",4]
mc[1,"exec.gov.coef"]<-summary(exec)$coefficients["gov.rep",1]
mc[1,"exec.gov.p"]<-summary(exec)$coefficients["gov.rep",4]

return(mc)

}

##########################################
#Write a function to create a first-difference data set
fd<-function(data){
  holder<-matrix(NA,ncol=ncol(data),nrow=length(state)*8) #fd data set will have 8 periods
  colnames(holder)<-colnames(data)
  for(i in 2:9){ #start with second period
    for(j in 0:48){#49 states
      holder[(i-1)+8*j,1]<-data[i+9*j,1]-data[(i-1)+9*j,1] #first 8 columns should be differenced
      holder[(i-1)+8*j,2]<-data[i+9*j,2]-data[(i-1)+9*j,2]
      holder[(i-1)+8*j,3]<-data[i+9*j,3]-data[(i-1)+9*j,3]
      holder[(i-1)+8*j,4]<-data[i+9*j,4]-data[(i-1)+9*j,4]
      holder[(i-1)+8*j,5]<-data[i+9*j,5]-data[(i-1)+9*j,5]
      holder[(i-1)+8*j,6]<-data[i+9*j,6]-data[(i-1)+9*j,6]
      holder[(i-1)+8*j,7]<-data[i+9*j,7]-data[(i-1)+9*j,7]
      holder[(i-1)+8*j,8]<-data[i+9*j,8]-data[(i-1)+9*j,8]
      holder[(i-1)+8*j,9]<-data[i+9*j,9] #copy gubernatorial party indicator
      holder[(i-1)+8*j,10]<-data[i+9*j,10] #copy session number
      holder[(i-1)+8*j,11]<-data[i+9*j,11] #copy state
    }
  }
  
  holder<-data.frame(holder, stringsAsFactors = F)
  
  for(t in 1:(ncol(holder)-1)){#convert to numeric
    holder[,t]<-as.numeric(holder[,t])
  }
  
  return(holder)
}

##########################################
#Write a function to run regressions and store coefficients and p-values

fd.reg<-function(data){
  
  names<-c("pivot.med.coef","pivot.med.p","pivot.maj.coef","pivot.maj.p","pivot.gov.coef","pivot.gov.p",
           "pivot.sm.med.coef","pivot.sm.med.p","pivot.sm.maj.coef","pivot.sm.maj.p","pivot.sm.gov.coef","pivot.sm.gov.p",
           "median.med.coef","median.med.p","median.maj.coef","median.maj.p","median.gov.coef","median.gov.p",
           "gatekeep.med.coef","gatekeep.med.p","gatekeep.maj.coef","gatekeep.maj.p","gatekeep.gov.coef","gatekeep.gov.p",
           "positive.med.coef","positive.med.p","positive.maj.coef","positive.maj.p","positive.gov.coef","positive.gov.p",
           "exec.med.coef","exec.med.p","exec.maj.coef","exec.maj.p","exec.gov.coef","exec.gov.p")
  
  #create a dataframe for the store regression coefficients
  mc<-matrix(NA,nrow=1,ncol=length(names))
  colnames(mc)<-names
  
  #estimate models
  pivot<-lm(p.pivot~xmedian+xmaj+gov.rep+as.factor(st)-1,data=data)
  pivot.small<-lm(p.pivot.small~xmedian+xmaj+gov.rep+as.factor(st)-1,data=data)
  median<-lm(p.median~xmedian+xmaj+gov.rep+as.factor(st)-1,data=data)
  gatekeep<-lm(p.gatekeep~xmedian+xmaj+gov.rep+as.factor(st)-1,data=data)
  positive<-lm(p.positive~xmedian+xmaj+gov.rep+as.factor(st)-1,data=data)
  exec<-lm(p.exec~xmedian+xmaj+gov.rep+as.factor(st)-1,data=data)
  
  #store coefficients and p-values
  mc[1,"pivot.med.coef"]<-summary(pivot)$coefficients["xmedian",1]
  mc[1,"pivot.med.p"]<-summary(pivot)$coefficients["xmedian",4]
  mc[1,"pivot.maj.coef"]<-summary(pivot)$coefficients["xmaj",1]
  mc[1,"pivot.maj.p"]<-summary(pivot)$coefficients["xmaj",4]
  mc[1,"pivot.gov.coef"]<-summary(pivot)$coefficients["gov.rep",1]
  mc[1,"pivot.gov.p"]<-summary(pivot)$coefficients["gov.rep",4]
  
  
  mc[1,"pivot.sm.med.coef"]<-summary(pivot.small)$coefficients["xmedian",1]
  mc[1,"pivot.sm.med.p"]<-summary(pivot.small)$coefficients["xmedian",4]
  mc[1,"pivot.sm.maj.coef"]<-summary(pivot.small)$coefficients["xmaj",1]
  mc[1,"pivot.sm.maj.p"]<-summary(pivot.small)$coefficients["xmaj",4]
  mc[1,"pivot.sm.gov.coef"]<-summary(pivot.small)$coefficients["gov.rep",1]
  mc[1,"pivot.sm.gov.p"]<-summary(pivot.small)$coefficients["gov.rep",4]
  
  mc[1,"median.med.coef"]<-summary(median)$coefficients["xmedian",1]
  mc[1,"median.med.p"]<-summary(median)$coefficients["xmedian",4]
  mc[1,"median.maj.coef"]<-summary(median)$coefficients["xmaj",1]
  mc[1,"median.maj.p"]<-summary(median)$coefficients["xmaj",4]
  mc[1,"median.gov.coef"]<-summary(median)$coefficients["gov.rep",1]
  mc[1,"median.gov.p"]<-summary(median)$coefficients["gov.rep",4]
  
  mc[1,"gatekeep.med.coef"]<-summary(gatekeep)$coefficients["xmedian",1]
  mc[1,"gatekeep.med.p"]<-summary(gatekeep)$coefficients["xmedian",4]
  mc[1,"gatekeep.maj.coef"]<-summary(gatekeep)$coefficients["xmaj",1]
  mc[1,"gatekeep.maj.p"]<-summary(gatekeep)$coefficients["xmaj",4]
  mc[1,"gatekeep.gov.coef"]<-summary(gatekeep)$coefficients["gov.rep",1]
  mc[1,"gatekeep.gov.p"]<-summary(gatekeep)$coefficients["gov.rep",4]
  
  mc[1,"positive.med.coef"]<-summary(positive)$coefficients["xmedian",1]
  mc[1,"positive.med.p"]<-summary(positive)$coefficients["xmedian",4]
  mc[1,"positive.maj.coef"]<-summary(positive)$coefficients["xmaj",1]
  mc[1,"positive.maj.p"]<-summary(positive)$coefficients["xmaj",4]
  mc[1,"positive.gov.coef"]<-summary(positive)$coefficients["gov.rep",1]
  mc[1,"positive.gov.p"]<-summary(positive)$coefficients["gov.rep",4]
  
  mc[1,"exec.med.coef"]<-summary(exec)$coefficients["xmedian",1]
  mc[1,"exec.med.p"]<-summary(exec)$coefficients["xmedian",4]
  mc[1,"exec.maj.coef"]<-summary(exec)$coefficients["xmaj",1]
  mc[1,"exec.maj.p"]<-summary(exec)$coefficients["xmaj",4]
  mc[1,"exec.gov.coef"]<-summary(exec)$coefficients["gov.rep",1]
  mc[1,"exec.gov.p"]<-summary(exec)$coefficients["gov.rep",4]
  
  return(mc)
  
}

########################################################

#Monte Carlo Simulation ###################################

b<-10000 #set the number of simulations

state<-unique(MCData$state) #get vector of state names for the sim function

# dopar setup for parallel computing
detectCores()
cl<-makePSOCKcluster(6)
registerDoParallel(cl)
clusterExport(cl, varlist=ls())

#note this will produce errors for each call of summary of the median model regression - should be 6xb;

#Get max and min for range of status quos
min(MCData$hou_dem.avg,na.rm=F)
max(MCData$hou_rep.avg,na.rm=F)

mc<-foreach(i = 1:b, .combine="rbind") %dopar% {
  sim.data<-sim(state=state)
  lv<-level.reg(data=sim.data)
  lv1<-matrix(NA,ncol=ncol(lv)+1,nrow=1)
  colnames(lv1)<-c(colnames(lv),"reg")
  lv1[,1:ncol(lv)]<-lv
  lv1[,ncol(lv)+1]<-"lv"
  
  fd.data<-fd(data=sim.data)
  fdf<-fd.reg(data=fd.data)
  fdf1<-matrix(NA,ncol=ncol(fdf)+1,nrow=1)
  fdf1[,1:ncol(fdf)]<-fdf
  fdf1[,ncol(fdf)+1]<-"fd"
  
  rbind(lv1,fdf1)
}

stopCluster(cl)

##########################
#create table of level reg results
##########################

lv.only<-mc[which(mc[,ncol(mc)]=="lv"),]
lv.only<-lv.only[,-ncol(mc)] #remove reg indicator

out<-matrix(NA,ncol=ncol(lv.only),nrow=1)
colnames(out)<-colnames(lv.only)

for(i in seq(1,(ncol(lv.only)-1),by=2)){ #get means of coefficients
  out[1,i]<-mean(as.numeric(lv.only[,i]))
}

for(i in seq(2,ncol(lv.only),by=2)){ #get proportion of coefficients that are positive and have a p-value <0.05
  x<-ifelse(as.numeric(lv.only[,i])<0.05 & as.numeric(lv.only[,(i-1)])>0,1,0)
  out[1,i]<-mean(x)
}

table.out.lvl<-rbind(out[1:6],out[7:12],out[13:18],out[19:24],out[25:30],out[31:36])

colnames(table.out.lvl)<-c("Chamber Median Avg Coef","prop. of chamber p-values < 0.05",
                       "Majority Median Avg Coef","prop. of majority p-values < 0.05",
                       "Rep. Gov Indicator Avg Coef","prop. of governor p-values < 0.05")

rownames(table.out.lvl)<-c("Pivot","Pivot No Filibuster", "Median", "Negative Agenda Power", "Positive Agenda Power", "Exec. Agenda Setter")


##########################
#create table of fd reg results
##########################

fd.only<-mc[which(mc[,ncol(mc)]=="fd"),]
fd.only<-fd.only[,-ncol(mc)] #remove reg indicator

out<-matrix(NA,ncol=ncol(fd.only),nrow=1)
colnames(out)<-colnames(fd.only)

for(i in seq(1,(ncol(fd.only)-1),by=2)){ #get means of coefficients
  out[1,i]<-mean(as.numeric(fd.only[,i]))
}

for(i in seq(2,ncol(fd.only),by=2)){ #get proportion of coefficients that have a p-value <0.05
  x<-ifelse(as.numeric(fd.only[,i])<0.05 & as.numeric(fd.only[,(i-1)])>0,1,0)
  out[1,i]<-mean(x)
}

table.out.fd<-rbind(out[1:6],out[7:12],out[13:18],out[19:24],out[25:30],out[31:36])

colnames(table.out.fd)<-c("Chamber Median Avg Coef","prop. of chamber p-values < 0.05",
                       "Majority Median Avg Coef","prop. of majority p-values < 0.05",
                       "Rep. Gov Indicator Avg Coef","prop. of governor p-values < 0.05")

rownames(table.out.fd)<-c("Pivot","Pivot No Filibuster", "Median", "Negative Agenda Power", "Positive Agenda Power", "Exec. Agenda Setter")



##########################
#Assume all states require a 3/5 supermajority to invoke cloture
##########################

#set fp = 3/5 in sim() function

sim35<-function(state){
  
  #create a matrix for the simulated data
  data<-matrix(NA,nrow=1,ncol=11)
  colnames(data)<-c("xmaj","xmedian","p.pivot","p.pivot.small","p.median","p.gatekeep","p.positive","p.exec","gov.rep","session","st")
  
  #Loop to model each state session and create a simulated data set
  for(j in 1:length(state)){
    
    sq<-runif(1,min=-3,max=3)
    
    #Get data for the state being modeled
    dem.avg<-MC35$hou_dem.avg[MC35$state==state[j]]
    dem.sd<-MC35$hou_dem.sd[MC35$state==state[j]]
    rep.avg<-MC35$hou_rep.avg[MC35$state==state[j]]
    rep.sd<-MC35$hou_rep.sd[MC35$state==state[j]]
    length.rep<-MC35$rep.size[MC35$state==state[j]]
    gov.rep<-MC35$gov.rep[MC35$state==state[j]]
    vp<-unique(MC35$vp[MC35$state==state[j]])
    fp<-3/5
    sessions<-length(dem.avg)
    
    #get the simulated sessions for the state
    x<-lawmaking(dem.avg=dem.avg,dem.sd=dem.sd,rep.avg=rep.avg,rep.sd=rep.sd,gov.rep=gov.rep,length.rep=length.rep,vp=vp,fp=fp,sq=sq,sessions=sessions)
    
    st<-rep(state[j],length=sessions)#create a state label
    
    session<-MC35$session[MC35$state==state[j]] #get session number for regressions
    
    x<-cbind(x,session,st) #combine the state data, state label, and session number 
    
    data<-rbind(data,x) #add the state to the full data set
  }
  
  data<-as.data.frame(data,stringsAsFactors=F) #convert to a data frame
  
  #convert variables to numeric
  for(t in 1:(ncol(data)-1)){
    data[,t]<-as.numeric(data[,t])
  }
  data<-data[-1,] #remove the first row of NA's used to start the data 
  return(data)
}

# dopar setup for parallel computing
detectCores()
cl<-makePSOCKcluster(6)
registerDoParallel(cl)
clusterExport(cl, varlist=ls())

#run the simulation
mc35<-foreach(i = 1:b, .combine="rbind") %dopar% {
  sim.data<-sim35(state=state)
  lv<-level.reg(data=sim.data)
  lv1<-matrix(NA,ncol=ncol(lv)+1,nrow=1)
  colnames(lv1)<-c(colnames(lv),"reg")
  lv1[,1:ncol(lv)]<-lv
  lv1[,ncol(lv)+1]<-"lv"
  
  fd.data<-fd(data=sim.data)
  fdf<-fd.reg(data=fd.data)
  fdf1<-matrix(NA,ncol=ncol(fdf)+1,nrow=1)
  fdf1[,1:ncol(fdf)]<-fdf
  fdf1[,ncol(fdf)+1]<-"fd"
  
  rbind(lv1,fdf1)
}

stopCluster(cl)

#get table of level results for 3/5 simulation
lv.35<-mc35[which(mc35[,ncol(mc35)]=="lv"),]
lv.35<-lv.35[,-ncol(mc35)] #remove reg indicator

out<-matrix(NA,ncol=ncol(lv.35),nrow=1)
colnames(out)<-colnames(lv.35)

for(i in seq(1,(ncol(lv.35)-1),by=2)){ #get means of coefficients
  out[1,i]<-mean(as.numeric(lv.35[,i]))
}

for(i in seq(2,ncol(lv.35),by=2)){ #get proportion of coefficients that are positive and have a p-value <0.05
  x<-ifelse(as.numeric(lv.35[,i])<0.05 & as.numeric(lv.35[,(i-1)])>0,1,0)
  out[1,i]<-mean(x)
}

rw<-rownames(table.out.lvl)
table.out.lvl<-rbind(table.out.lvl,out[1:6])
rownames(table.out.lvl)<-c(rw,"Pivot Model (All are 3/5)")

#create table of fd results for 3/5 simulation
fd.35<-mc35[which(mc35[,ncol(mc35)]=="fd"),]
fd.35<-fd.35[,-ncol(mc35)] #remove reg indicator

out<-matrix(NA,ncol=ncol(fd.35),nrow=1)
colnames(out)<-colnames(fd.35)

for(i in seq(1,(ncol(fd.35)-1),by=2)){ #get means of coefficients
  out[1,i]<-mean(as.numeric(fd.35[,i]))
}

for(i in seq(2,ncol(fd.35),by=2)){ #get proportion of coefficients that have a p-value <0.05
  x<-ifelse(as.numeric(fd.35[,i])<0.05 & as.numeric(fd.35[,(i-1)])>0,1,0)
  out[1,i]<-mean(x)
}

rw<-rownames(table.out.fd)
table.out.fd<-rbind(table.out.fd,out[1:6])
rownames(table.out.fd)<-c(rw,"Pivot Model (All are 3/5)")

##########################
#Assume certain states require unanimity
##########################

#set fp = MCunam$fp in sim() function

sim_unam<-function(state){
  
  #create a matrix for the simulated data
  data<-matrix(NA,nrow=1,ncol=11)
  colnames(data)<-c("xmaj","xmedian","p.pivot","p.pivot.small","p.median","p.gatekeep","p.positive","p.exec","gov.rep","session","st")
  
  #Loop to model each state session and create a simulated data set
  for(j in 1:length(state)){
    
    sq<-runif(1,min=-3,max=3)
    
    #Get data for the state being modeled
    dem.avg<-MCunam$hou_dem.avg[MCunam$state==state[j]]
    dem.sd<-MCunam$hou_dem.sd[MCunam$state==state[j]]
    rep.avg<-MCunam$hou_rep.avg[MCunam$state==state[j]]
    rep.sd<-MCunam$hou_rep.sd[MCunam$state==state[j]]
    length.rep<-MCunam$rep.size[MCunam$state==state[j]]
    gov.rep<-MCunam$gov.rep[MCunam$state==state[j]]
    vp<-unique(MCunam$vp[MCunam$state==state[j]])
    fp<-unique(MCunam$fp[MCunam$state==state[j]])
    sessions<-length(dem.avg)
    
    #get the simulated sessions for the state
    x<-lawmaking(dem.avg=dem.avg,dem.sd=dem.sd,rep.avg=rep.avg,rep.sd=rep.sd,gov.rep=gov.rep,length.rep=length.rep,vp=vp,fp=fp,sq=sq,sessions=sessions)
    
    st<-rep(state[j],length=sessions)#create a state label
    
    session<-MCunam$session[MCunam$state==state[j]] #get session number for regressions
    
    x<-cbind(x,session,st) #combine the state data, state label, and session number 
    
    data<-rbind(data,x) #add the state to the full data set
  }
  
  data<-as.data.frame(data,stringsAsFactors=F) #convert to a data frame
  
  #convert variables to numeric
  for(t in 1:(ncol(data)-1)){
    data[,t]<-as.numeric(data[,t])
  }
  data<-data[-1,] #remove the first row of NA's used to start the data 
  return(data)
}


# dopar setup for parallel computing
detectCores()
cl<-makePSOCKcluster(6)
registerDoParallel(cl)
clusterExport(cl, varlist=ls())

#run the simulation
mc_unam<-foreach(i = 1:b, .combine="rbind") %dopar% {
  sim.data<-sim_unam(state=state)
  lv<-level.reg(data=sim.data)
  lv1<-matrix(NA,ncol=ncol(lv)+1,nrow=1)
  colnames(lv1)<-c(colnames(lv),"reg")
  lv1[,1:ncol(lv)]<-lv
  lv1[,ncol(lv)+1]<-"lv"
  
  fd.data<-fd(data=sim.data)
  fdf<-fd.reg(data=fd.data)
  fdf1<-matrix(NA,ncol=ncol(fdf)+1,nrow=1)
  fdf1[,1:ncol(fdf)]<-fdf
  fdf1[,ncol(fdf)+1]<-"fd"
  
  rbind(lv1,fdf1)
}

stopCluster(cl)

#get table of level results for unanimity
lv.unam<-mc_unam[which(mc_unam[,ncol(mc_unam)]=="lv"),]
lv.unam<-lv.unam[,-ncol(mc_unam)] #remove reg indicator

out<-matrix(NA,ncol=ncol(lv.unam),nrow=1)
colnames(out)<-colnames(lv.unam)

for(i in seq(1,(ncol(lv.unam)-1),by=2)){ #get means of coefficients
  out[1,i]<-mean(as.numeric(lv.unam[,i]))
}

for(i in seq(2,ncol(lv.unam),by=2)){ #get proportion of coefficients that are positive and have a p-value <0.05
  x<-ifelse(as.numeric(lv.unam[,i])<0.05 & as.numeric(lv.unam[,(i-1)])>0,1,0)
  out[1,i]<-mean(x)
}

rw<-rownames(table.out.lvl)
table.out.lvl<-rbind(table.out.lvl,out[1:6])
rownames(table.out.lvl)<-c(rw,"Pivot Model (Unanimity)")

#create table of fd results for unanimity
fd.unam<-mc_unam[which(mc_unam[,ncol(mc_unam)]=="fd"),]
fd.unam<-fd.unam[,-ncol(mc_unam)] #remove reg indicator

out<-matrix(NA,ncol=ncol(fd.unam),nrow=1)
colnames(out)<-colnames(fd.unam)

for(i in seq(1,(ncol(fd.unam)-1),by=2)){ #get means of coefficients
  out[1,i]<-mean(as.numeric(fd.unam[,i]))
}

for(i in seq(2,ncol(fd.unam),by=2)){ #get proportion of coefficients that have a p-value <0.05
  x<-ifelse(as.numeric(fd.unam[,i])<0.05 & as.numeric(fd.unam[,(i-1)])>0,1,0)
  out[1,i]<-mean(x)
}

rw<-rownames(table.out.fd)
table.out.fd<-rbind(table.out.fd,out[1:6])
rownames(table.out.fd)<-c(rw,"Pivot Model (Unanimity)")

#############################################################
#save tables of results and sumlation data
#############################################################

save(mc,file="MCregs.RData")
save(mc35,file="MCregs_three_fifths.RData")
save(mc_unam,file="MCregs_unanimity.RData")

write.csv(round(table.out.lvl,digits=4),file="MCsummary_level.csv")
write.csv(round(table.out.fd,digits=4),file="MCsummary_fd.csv")
